package it.uniroma2.dtk.op.convolution;

import java.util.Arrays;

import it.uniroma2.dtk.op.TransformAndCompose;
import it.uniroma2.util.math.Complex;
import it.uniroma2.util.math.FFT;
import it.uniroma2.util.vector.VectorProvider;

public class CircularConvolution extends TransformAndCompose {

	protected FFT fft = null;		// Object used to optionally compute the Fast Fourier Transform for the circular convolutions
	protected int spaceGrade = 0;	// log_2 of the vector dimension, where an integer, for using the FFT
	
	@Override
	public double[] op(double[] x, double[] y) {
		if (x == null)
			return y;
		else if (y == null)
			return x;
		if (fft == null)
			return circularConvolutionBasic(x, y);
		else
			return circularConvolutionFFT(x, y);
	}
	
	protected double[] circularConvolutionBasic(double[] firstVector, double[] secondVector) {
		int size = firstVector.length;
		double[] result = new double[size];
		Arrays.fill(result, 0);
		for (int i=0; i<size; i++)
			for (int j=0; j<size; j++)
				result[i] += firstVector[j] * secondVector[(i-j)<0 ? (i-j+size) : i-j];
		return result;
	}
	
	protected double[] circularConvolutionFFT(double[] firstVector, double[] secondVector) {
	return Complex.extract(fft.real_cconvolve(
					Complex.generate(firstVector),
					Complex.generate(secondVector),
					spaceGrade));
	}
	
	@Override
	public void initialize(VectorProvider vp) throws Exception {
		super.initialize(vp);
		int exp = log2(vp.getVectorSize());
		if (exp < 0)
			return;
		else {
			fft = new FFT();
			fft.initialize(exp);
			spaceGrade = exp;
		}
	}
	
	protected int log2(int arg) {
		int exp = 0;
		while(true) {
			double pow = Math.pow(2, exp); 
			if (pow == arg)
				return exp;
			else if (pow > arg)
				return -1;
			exp++;
		}
	}

}
